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1.  Introduction 


The  U.S.  Army  has  a  growing  interest  in  the  use  of  advanced  sensors  and  computer  models  to 
retrieve,  display,  and  interpret  acoustic  signals  on  future  battlefields.  Such  technologies  are 
emerging  for  the  surveillance,  detection,  identification,  and  tracking  of  sound-emitting  targets  in 
and  around  forests,  hilly  terrain,  and  in  cities  (1).  Consequently,  the  U.S.  Army  is  looking  to 
implement  the  best  possible  computer  models  for  determining  point-to-point  acoustic 
transmission.  At  the  same  time,  the  retrieval  and  interpretation  of  acoustic  signals  are  greatly 
influenced  by  turbulence  and  refraction  effects  caused  by  finer  scale  atmospheric  motions  over 
varying  topography  and  surface  energy  budgets  (2).  Hence,  improved  physics-based  theory  and 
computer  models  for  meteorology  coupled  to  acoustics  are  expected  to  contribute  important 
information  on  the  performance  of  future  combat  (acoustic)  systems. 

Recently,  several  physics-based  computer  models  have  been  developed  to  calculate  one-  and 
two-dimensional  forest  canopy  micrometeorology  and  turbulence  for  future  U.S.  Army  acoustic 
application  research  (3,  4).  Individual  computer  codes  have  incorporated  various  computational 
schemes  on  a  uniform  grid  to  solve  the  meteorological  fields  (5).  However,  it  may  be  possible  to 
improve  the  efficiency  and  realism  of  the  coupled  meteorological — acoustic  computer  models  by 
introducing  variable  grid.  Variable  grid  will  allow  for  better  distribution  of  grid  points  and  will 
extend  calculations  higher  into  the  boundary  layer  above  the  forest.  In  addition,  a  finer  grid 
inside  the  forest  and  a  coarser  grid  above  the  forest  will  help  to  resolve  important  meteorological 
(and  acoustic)  scales  and  processes. 

Thus,  the  following  report  presents  results  from  some  preliminary  tests  to  incorporate  this 
important  feature  into  numerical  codes.  First,  several  simpler  diffusion  models  are  developed  to 
benchmark  fundamental  (numerical)  techniques.  Both  explicit  and  implicit  differencing  schemes 
are  examined.  In  addition,  numerical  stability  criteria  for  these  calculations  are  demonstrated. 
Then,  successful  preliminary  tests  on  these  codes  are  extended  to  more  complicated 
meteorological — acoustic  models  for  the  forest  canopy. 


2.  Explicit  Differencing 


A  simplified  conservation  equation  to  describe  the  mean  concentration  (diffusion)  of  a  scalar  C 
can  be  written  as  follows: 

dC  _  dC  d^' 

—  =  -  u - ,  ( 

dt  dx  dz 
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where  t  is  the  independent  variable  time,  u  is  the  mean  longitudinal  component  of  the  wind 
velocity,  x  is  range,  z  is  height  above  ground,  and  w'C  is  the  mean  scalar  flux.  The  flux- 
gradient  assumption  (6)  suggests  that 


-WC  =  K^  ,  (2) 

dz 


where  K  is  the  scalar  (eddy)  diffusivity.  Combining  equations  1  and  2  yields. 


dt 


—  —  u - h  K 

dx 


d^c 

dz^ 


(3) 


In  discretized  form,  this  simple  model  can  be  solved  forward  in  time  using  the  following  explicit 
differencing  scheme  for  variable  grid,  i.e.. 


=c"j  +At| 


^  n  ^  n 

■Hj— - -  +  ^- 


X,  -  X 


(z,  -  z,.,  ) + -  z,  >  -  (z,,,  -  z,., ) 

(z  M  -  Z  )(z  M  -  z  ,  )(z  .  -  Z  .  1 ) 


,  (4) 


where  i  and  j  are  the  indices  for  the  horizontal  and  vertical  grid,  respectively,  and  n  is  the  time- 
step.  To  execute  the  model,  the  following  steps  may  be  taken: 

1 .  Initialize  the  scalar  field  as  C] j  =  e  ,  where  h  defines  the  center  and  cr  defines  the 
width  of  the  Gaussian  profile. 

2.  Assume  the  wind  flow  is  constant,  i.e.,  u^  j  =3.0  . 

3.  Let  the  horizontal  grid  be  constant,  i.e..  Ax  =  0.5  m. 

4.  Assume  horizontal  scalar  gradients  are  zero,  i.e.,  apply  the  explicit  boundary  condition 
C2J  =  C"j  at  the  upstream  lateral  edge. 

5.  Implement  the  following  variable  grid: 

Az  =  0.25  m  0.25  <z<  10.0  m  , 

Az  =  0.50  m  10.5  <z<  30.0  m  , 

Az=  1.00m  31.0<z<  50.0m  . 

Figure  1  shows  the  model  results  for  the  case  where  the  scalar  diffusivity  is  A  =  10  m  s  and  the 
time-step  is  At  =  0.001  s.  These  are  the  time-height  contours  derived  from  a  simple  scalar 
diffusion  model  solved  via  an  explicit  differencing  scheme.  Here,  the  time-step  is  small  to 
satisfy  the  stability  criterion  described  by  Press  et  al.  (7)  as 
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Figure  1.  Time -height  contours  derived  from  a  scalar  diffusion  model 
solved  via  an  explicit  differencing  scheme.  In  this  example, 
the  scalar  diffusivity  is  10  s' '  and  the  time-step  is 
At  =  0.001  s. 


2KAt 

w 


<1 


(5) 


Otherwise,  for  larger  time-steps,  the  numerieal  seheme  would  be  unstable  and  the  model  would 
not  produee  viable  results.  Alternately,  figure  2  shows  several  profiles  of  sealar  eoneentration  at 
different  time-steps  derived  from  this  model.  Finally,  for  eomparison,  figure  3  shows  the  time- 
height  eontours  for  the  ease  where  K=3  m  s  and  the  time-step  is  Af  =  0.0 1  s. 


3.  Time-Dependent  Implicit  Differencing 


Alternately,  equation  3  (without  adveetion)  ean  be  solved  via  implieit  differeneing.  The 
diseretized  form  for  equation  3  in  this  ease  ean  be  written  as 


=KAt 


hj 


(6) 


3 


Time  (s) 


Figure  3.  Same  as  figure  1  except  the  scalar  diffusivity  K^3  s  * 
and  the  time-step  is  At  =  0.01  s. 
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If  it  is  assumed  that 


a  = 


_  KM  _ 

X  hj  -  )  ’ 


(V) 


then  equation  6  ean  be  solved  for  C.".  as 


-  a(z.  -  z  +  [i + «(z  .^1  -  z._, )]  c,;;'  -  «(z 


-^7 


n+\  _  ^  n 

*-'1,7-1  “  ‘-'1,7 


(8) 


Now,  equation  8  can  be  solved  via  a  tridiagonal  matrix  algorithm  (also  called  the  Thomas 
algorithm)  (Press  et  al.  [5]).  By  implementing  the  Thomas  algorithm,  one  is  looking  for  a 
solution  to  the  equation  A  •  x  =  r,  which  can  be  written  in  matrix-vector  notation  as 


1 

o 

1 

1 

1 _ _ 

1 

1 _ — 

a2  hj  Cj 

c^_[ 

^N-\ 

: 

1 

0  J 

1 

1 _ 

1 

_ 1 

(9) 


Here  a,  b,  c,  and  r  are  the  input  vectors  and  x  is  the  (scalar)  output  vector.  The  sub-,  main-,  and 
super-diagonal  vectors  (a,  b,  and  c,  respectively)  can  be  computed  from  equation  8  as 

a,j=-cc[zj^x-Zj)  ,  (10) 

=l  +  a(z.^i-z._j)  ,  (11) 

and 

•  (12) 


The  input  vector  r  is  defined  as. 


KJ 


=  C, 


^7 


(13) 


In  addition,  the  following  boundary  conditions  are  implemented:  atzi,  a/j  =  0,  bi_i  =  1,  Ci_\  =  0, 
and  r.  y  =  C,.  j  and  at  z^,  ai,Ar=  0,  bi_N=  1,  c,'7v=  0,  and  r.  f^  =  C,  ^ . 


Figure  4  shows  the  time-height  contours  derived  from  this  scalar  diffusion  model  solved  via 
implicit  differencing,  where  the  scalar  diffusivity  is  A  =  10  m  s  and  the  time-step  is  At  =  0. 1  s. 
As  expected,  these  results  are  quite  similar  to  those  shown  earlier  derived  via  explicit 
differencing  (figure  1).  Except  here,  the  magnitude  of  the  time-step  is  not  critical  because  the 
stability  criterion  for  this  model  is  given  instead  as 


5 


(14) 


1  + 


(Azf 


<1 


so  that  the  numerical  scheme  is  unconditionally  stable  for  both  large  and  small  At  (P). 

Alternately,  figures  5  and  6  show,  for  comparison,  the  time -height  contours  for  the  case  where  At 
=  0.5  5  and  At  =  1.0  5,  respectively. 


Figure  4.  Time-height  contours  derived  from  a  scalar  diffusion  model 
solved  via  an  implicit  differencing  scheme.  In  this  example, 
the  scalar  diffusivity  is  =  10  ^  *  and  the  time-step  is 

At  =  0.1 
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Figure  6.  Same  as  figure  4  except  the  time-step  is  At  =  1.0  s. 


1 


4.  Steady-State  Implicit  Differencing 


4,1  Scalar  Diffusion  Model 

In  contrast,  a  steady-state  model  ean  be  derived  from  equation  3,  i.e., 

d^C_  u  dC 
dz^  K  dx 

whieh  can  be  written  in  diseretized  form  as, 

(^y  -  )+  ^ 


If  it  is  assumed  that 


and  /3  =  y2  (z  .,1  -  z  .^1  “  ^y-i ) 


then  equation  16  ean  be  solved  for  C,_i  .  as 


(z  .  -  Z  )  c,.  +  1  +  ^  (^yU  -  ^yU  )  C’-.y  “  ^  tyU  “  ^y  )  <^.-,y-l  =  C’/-!,./  ' 


Here,  equation  18  ean  be  solved  via  the  Thomas  algorithm,  where  the  input  vectors  a,  b,  c,  and  r 


are  given  as 


a  i 


^/,y  =l  +  ^tyU-Vi) 

^/,y  =-^(^y  -Vi)  ’ 


'i.y  =  ■  (2: 

Figure  7  shows  the  two-dimensional  (2-D)  sealar  field  as  eomputed  by  the  steady-state  model. 
For  this  example,  the  vertieal  and  horizontal  grid  dimensions  (i.e.,  Az  =  0.25,  0.5,  1.0  m  and  Ax  ■ 
0.5  m)  and  the  initial  sealar  profile  are  the  same  as  those  deseribed  earlier  in  the  report.  In 
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contrast,  figure  8  shows  the  2-D  steady-state  sealar  field  for  the  case  where  Ax  =  2.0  m.  Here, 
the  magnitude  of  the  horizontal  grid  increment  is  not  eritical  because  the  stability  criterion  for 
this  model  is  given  as 


1  + 


AK^t 

(Azf 


<1 


(23) 


where  one  may  substitute  At  =  Ax/m  .  Thus,  a  simple,  unconditionally  stable  2-D  diffusion 
model  has  been  developed  to  test  the  basic  formulation  and  correct  implementation  of  variable 
grid  and  implicit  differencing  to  solve  the  steady  state  sealar  field. 


Figure  7.  The  2-D  steady-state  scalar  field  derived  from  a  scalar 

diffusion  model  solved  via  an  implicit  differencing  scheme. 

In  this  example,  the  scalar  diffusivity  ^  '  and 

the  horizontal  grid  is  Ax  =  0.5  m. 

4.2  Coupled  Meteorological — ^Acoustic  Computer  Model 

In  this  section,  successful  preliminary  tests  on  simpler  (diffusion)  codes  are  extended  to  more 
complicated  meteorological — acoustic  models  for  the  forest  canopy.  As  an  example,  Tunick  (4) 
gives  the  parameterized  equation  for  Reynolds  stress  (u’w')  as 


9 


where  is  the  mean  flow  -  longitudinal,  (wj  is  the  mean  flow  -  vertieal,  (u'^j  and  are 

the  longitudinal  and  vertieal  veloeity  varianees,  respeetively,  (u'0'^  is  the  horizontal  heat  flux,  g 

is  the  aeeeleration  due  to  gravity,  T  is  the  absolute  temperature,  C  is  a  eonstant  whose  value  is 
~0.077,  q  is  the  turbulent  kinetie  energy,  Cd  is  the  forest  eanopy  drag  eoeffieient,  A  (in  units 
m  m'  )  is  the  leaf  area  density,  and  A  is  a  funetion  of  mixing  length^  The  overbar  and  primed 

hn  equation  24,  /Ij  and  A2  are  length  scales,  which  contain  two  necessary  closure  constants,  i.e.,  1^-  =  ,  where  k  is  an 

arbitrary  index  and  I  is  the  mixing  length.  Values  for  these  closure  constants  are  =  0.39  and  02  =  0.85  {10). 
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variables  indicate  the  mean  (time  averaged)  and  fluctuating  components  of  the  given  quantity, 
whereas  the  brackets,  ^  ^  ,  indicate  horizontal  averaging.  Solving  equation  (24)  for  d(^  jdz 

yields 


—\\d{u 
w'  , , 

dz 


=  {u 


d(u'w')  d(u'w')  . _ ,  dlu)  / — \  s/w)  , _ >  d{w 

.A  \  /  ,  A.A  \  /  +  \  ' 


A 

dx 


dx 
f  f 

qX, 

V  V 


-  +  (w 


dz 


d(u'w'^  diu 


dx 


dz 


A 

dz 


J) 


dx 
f  f 

(iK 

V  V 


dx 


I  dx  ''  /  dz 


JJ 


d(w'^j  d(u'w'^ 


dz 


-^(ud')  +  -^(u'w' 

t\  /  3A. 


diw) 

-Cq^  ^^-{w)C„A 
dx 


U\ 


u)-{u)C^A 


U\ 


w 


(25) 


so  that 


dz 


=  GRADU, , 

GJ 


(26) 


where  r.h.s.  contains  all  the  terms  on  the  right-hand  side  of  equation  25.  Based  on  this 
expression,  the  following  second-order  ordinary  differential  equation  can  be  set  up  to  solve  for 


d^iu) 


^Mgradu,]  . 


dz  dz 

Equation  27  can  be  discretized  over  a  variable  grid  as 


(27) 


[z.  -[z.^,  +(. 

0-5 An  -  Vi  A'u  -^jhl  -  Vi) 


,  _  GRADU,^,,  - GRADU^^^_, 


.  (28) 


Now,  equation  28  can  be  solved  via  the  tridiagonal  matrix  algorithm,  where  the  sub-,  main-,  and 
super-diagonal  vectors  (a,  b,  and  c,  respectively)  are  as  follows: 
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= 


ij-l 


(29) 


-{u 

jj  -  '  '  ‘’J _ 

-^jhj  -  vi)  ’ 


^7 


(30) 


and 


c 


(31) 


Finally,  the  input  vector  r  is  given  as 

GRADU.  -  GRADU.  ,  , 

r,,.= - ^  .  (32) 

This  procedure  is  similarly  applied  to  solve  the  other  2-D  meteorological  fields,  which  include 
mean  temperature,  vertical  velocity,  fiuctuation  pressure,  and  heat  flux. 

Figure  9  shows  the  computed  horizontal  wind  velocity  ,  in  units  ms'\  over  a  variable  grid 

within  and  above  a  uniform  forest  stand.  For  this  example,  canopy  height  is  10  m.  In  contrast, 
figure  10,  shows  the  computed  horizontal  wind  velocity,  ,  and  wind  flow  streamlines  over  a 

variable  grid  within  and  above  a  non-uniform  forest  stand,  i.e.,  one  that  contains  multiple  step 
changes  in  canopy  height  at  the  lower  boundary. 


5.  Summary 


Introducing  variable  grid  will  improve  the  efficiency  and  realism  of  coupled  meteorological — 
acoustic  computer  models.  Variable  grid  will  allow  for  better  distribution  of  grid  points  and  will 
extend  calculations  higher  into  the  boundary  layer  above  the  forest.  A  finer  grid  inside  the  forest 
and  a  coarser  grid  above  the  forest  will  help  to  resolve  important  meteorological  (and  acoustic) 
scales  and  processes.  Therefore,  this  report  has  presented  results  from  some  preliminary  tests  to 
incorporate  variable  grid  into  current  software.  Several  basic  model  codes  were  developed  to 
benchmark  the  fundamental  (numerical)  techniques,  i.e.,  to  solve  diffusion  equations  via  explicit 
and  implicit  differencing.  These  calculations  worked  quite  well.  In  addition,  numerical  stability 
criteria  for  these  calculations  were  clearly  demonstrated. 
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Figure  9.  Fiorizontal  wind  velocity,  (u^  ,  in  units  ms  ',  within  and 

above  a  uniform  forest  stand,  where  canopy  height  Qi)  is 
10  m. 


Although  difficult  to  perform,  stability  analysis  is  important  so  that  models  will  be 
eomputationally  robust.  In  this  report,  it  was  shown  that  the  magnitude  of  the  time-step  (and/or 
horizontal  grid  inerement)  was  not  oritieal  for  the  time-dependent  and  steady-state  diffusion 
models  solved  via  implieh  differencing.  However,  in  eurrent  meteorologieal — acoustie  software, 
the  stable  ealeulation  of  the  seeond  horizontal  derivatives  (as  well  as  the  mixed  derivatives)  is 
found  to  neeesshate  fairly  large  Ax  (e.g.,  30-50  m).  Henee,  further  analysis  and  interpretation  of 
model  results  are  in  progress  to  resolve  some  remaining  difficulties  in  establishing  a  finer 
(horizontal)  grid  in  existing  meteorologieal — aeoustie  eodes. 
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Figure  10.  Florizontal  wind  velocity,  (u^  ,  in  units  ms'*,  and  the 

wind  flow  streamlines  within  and  above  a  non-uniform 
forest  stand,  i.e.,  one  that  contains  multiple  step  changes 
in  canopy  height  Qi)  at  the  lower  boundary. 
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